Prognosemodel - Uitval na 1 jaar - Zonder P

GVS | B Opleiding tot Verpleegkundige (HBO-V) - voltijd - versie 1.0

Auteur

Theo Bakker, lector Learning Technology & Analytics, De HHs

Publicatiedatum

6 juni 2024

1 Inleiding

1.1 Het nut van prognosemodellen

Prognosemodellen kunnen inzicht bieden in de factoren die gecorreleerd zijn aan de uitval van studenten. Met deze inzichten kan een opleiding interventies ontwikkelen om de uitval te verminderen of te voorkomen. Voor het lectoraat is het van belang dat we in vervolg hierop kunnen analyseren of er sprake is van bias in de data in relatie tot uitval en er mogelijk sprake is van een gebrek aan kansengelijkheid. Een prognosemodel is dus de opmaat naar een fairness analyse. Zie voor een verdere toelichting het onderzoeksprogramma ‘No Fairness without Awareness’.

1.2 Toelichting op de methode

Voor de ontwikkeling van prognosemodellen gebruiken we de aanpak van Tidymodels. Tidymodels is een framework voor het bouwen van een prognosemodel. Hiermee verzekeren we ons van een systematische, herhaalbare en schaalbare aanpak voor het bouwen van prognosemodellen.

1.3 Toelichting op de data

De basis voor deze analyse is studiedata van De Haagse Hogeschool (De HHs), verrijkt door het lectoraat Learning Technology & Analytics. De data bevat informatie over de inschrijvingen van studenten in het eerste jaar van de opleiding:

  1. Demografische kenmerken: geslacht, leeftijd, reistijd en SES totaalscore.
  2. Vooropleidingskenmerken: toelaatgevende vooropleiding, studiekeuzeprofiel, gemiddeld eindcijfer in de vooropleiding en eventuele deelname aan het Navitas programma.
  3. Aanmeldingskenmerken: aansluiting (direct na diploma, tussenjaar, switch), dag van aanmelding, aantal parallelle studies aan De HHs en collegejaar.

1.4 Toelichting op de analyse

We toetsen in deze analyse Uitval na 1 jaar bij studenten zonder propedeusediploma, voortaan Uitval genoemd.

Uitval is gedefinieerd als het niet meer ingeschreven staan in de opleiding in een aansluitend collegejaar. Een wisseling van opleidingsvorm binnen de opleiding, bijv. van voltijd in jaar 1 naar duaal in jaar 2, geldt niet als uitval.

2 Voorbereidingen

2.1 Laad de data

We laden een subset in van historische data specifiek voor:

Opleiding: GVS | B Opleiding tot Verpleegkundige (HBO-V), voltijd, eerstejaars - Uitval na 1 jaar

Toon code 
## Laad de data voor de HBO-V opleiding
dfOpleiding_inschrijvingen_base <- get_lta_studyprogram_enrollments_pin(
    board = "HHs/Inschrijvingen",
    faculty = faculteit,
    studyprogram = opleidingsnaam_huidig,
    studytrack = opleiding,
    studyform = toupper(opleidingsvorm),
    range = "eerstejaars")

## Herschik de levels
Set_Levels(dfOpleiding_inschrijvingen_base)

dfOpleiding_inschrijvingen_base <- dfOpleiding_inschrijvingen_base |>  
  
  ## Maak een eenvoudige Uitval variabele aan
  Mutate_Uitval(sUitval_model) |>
  
  ## Maak van de uitval variabele een factor
  mutate(SUC_Uitval = as.factor(SUC_Uitval)) |> 

  ## Verbijzonder eventueel op basis van het propedeusediploma
  Filter_Propedeusediploma(sPropedeusediploma) |>

  ## Maak van de Dubbele studie variabele een Ja/Nee variabele
  mutate(INS_Dubbele_studie = ifelse(INS_Aantal_inschrijvingen > 1, "Ja", "Nee")) |>  

  ## Verwijder INS_Aantal_inschrijvingen
  select(-INS_Aantal_inschrijvingen) |> 

  ## Pas voor een aantal variabelen de levels aan
  Mutate_Levels(
  c(
    "VOP_Studiekeuzeprofiel_LTA_afkorting",
    "INS_Aansluiting_LTA",
    "VOP_Toelaatgevende_vooropleiding_soort"
  ),
    list(lLevels_skp, lLevels_vop, lLevels_vop)
  )
  
## B Huidtherapie: Filter op uitsluitend studenten met een rangnummer (selectie)
if(opleiding == "HDT") {
  dfOpleiding_inschrijvingen_base <- dfOpleiding_inschrijvingen_base |> 
    filter(!is.na(RNK_Rangnummer)) 
}

2.2 Selecteer en inspecteer de data

We selecteren eerst de relevante variabelen. We verwijderen daarbij variabelen die maar 1 waarde hebben. We bekijken de variabelen in een samenvatting in relatie tot de Uitval. Daarnaast bekijken we de kwaliteit van de data op missende waarden.

Toon code 
lSelect <- c(
    "INS_Student_UUID_opleiding_vorm",
    "CBS_APCG_tf",
    "DEM_Geslacht",
    "DEM_Leeftijd_1_oktober",
    "GIS_Tijd_fiets_OV",
    "INS_Collegejaar",
    "INS_Dagen_tussen_aanmelding_en_1_september",
    "INS_Dubbele_studie",
    "INS_Aansluiting_LTA",
    "INS_Navitas_tf",
    "SES_Deelscore_arbeid",
    "SES_Deelscore_welvaart",
    "SES_Totaalscore",
    "SUC_Uitval",
    "VOP_Gemiddeld_cijfer_cijferlijst",
    "VOP_Gemiddeld_eindcijfer_VO_van_de_hoogste_vooropleiding_voor_het_HO",
    "VOP_Cijfer_CE1_nederlands",
    "VOP_Cijfer_CE1_engels",
    "VOP_Cijfer_CE_proxy_wiskunde",
    "VOP_Cijfer_CE1_natuurkunde",
    "VOP_Studiekeuzeprofiel_LTA_afkorting",
    "VOP_Toelaatgevende_vooropleiding_soort"
  )

## B Huidtherapie: voeg de variabele RNK_Rangnummer toe
if(opleiding == "HDT") {
  lSelect <- c(lSelect, "RNK_Rangnummer")
}

## Maak een subset
dfOpleiding_inschrijvingen <- dfOpleiding_inschrijvingen_base |>
  
  ## Selecteer de relevante variabelen
  select_at(lSelect) |>
  
  ## Hernoem variabelen voor beter leesbare namen
  rename(
    ID                    = INS_Student_UUID_opleiding_vorm,
    Geslacht              = DEM_Geslacht,
    Leeftijd              = DEM_Leeftijd_1_oktober,
    Reistijd              = GIS_Tijd_fiets_OV,
    Dubbele_studie        = INS_Dubbele_studie,
    Collegejaar           = INS_Collegejaar,
    Aanmelding            = INS_Dagen_tussen_aanmelding_en_1_september,
    Aansluiting           = INS_Aansluiting_LTA,
    Navitas               = INS_Navitas_tf,
    APCG                  = CBS_APCG_tf,
    SES_Arbeid            = SES_Deelscore_arbeid,
    SES_Welvaart          = SES_Deelscore_welvaart,
    SES_Totaal            = SES_Totaalscore,          
    Uitval                = SUC_Uitval,
    Cijfer_SE_VO          = VOP_Gemiddeld_cijfer_cijferlijst,
    Cijfer_CE_VO          = VOP_Gemiddeld_eindcijfer_VO_van_de_hoogste_vooropleiding_voor_het_HO,
    Cijfer_CE_Nederlands  = VOP_Cijfer_CE1_nederlands,
    Cijfer_CE_Engels      = VOP_Cijfer_CE1_engels,
    Cijfer_CE_Wiskunde    = VOP_Cijfer_CE_proxy_wiskunde,
    Cijfer_CE_Natuurkunde = VOP_Cijfer_CE1_natuurkunde,
    Studiekeuzeprofiel    = VOP_Studiekeuzeprofiel_LTA_afkorting,
    Vooropleiding         = VOP_Toelaatgevende_vooropleiding_soort
  ) |> 
  
  ## Pas CBS_APCG_tf aan naar factor
  mutate(APCG = case_when(APCG == TRUE ~ "Ja",
                          APCG == FALSE ~ "Nee",
                          .default = "Onbekend")) |>

  ## Geef aan waar missende cijfers in het VO zijn
  Mutate_Cijfers_VO() |>
  
  ## Verwijder variabelen, waarbij er maar 1 waarde is
  select(where(~ n_distinct(.) > 1)) |>
  
  arrange(Collegejaar, ID)

## B Huidtherapie: hernoem de variabele RNK_Rangnummer
if(opleiding == "HDT") {
  dfOpleiding_inschrijvingen <- dfOpleiding_inschrijvingen |> 
    rename(Rangnummer = RNK_Rangnummer)
} 

dfOpleiding_inschrijvingen <- dfOpleiding_inschrijvingen |> 
 ltabase::sort_distinct()

## Verwijder de basis dataset
rm(dfOpleiding_inschrijvingen_base)
Studentkenmerken versus Uitval
Variabele Uitval p-value2 Totaal, N = 18181
Ja, N=727 (40%)1 Nee, N=1091 (60%)1
Aanmelding 133,20 (74,18) 153,98 (72,88) <0,001 145,67 (74,08)
Aansluiting

0,025
    2e Studie 6 (32%) 13 (68%)
19 (100%)
    Direct 357 (38%) 579 (62%)
936 (100%)
    Na CD 4 (16%) 21 (84%)
25 (100%)
    Overig 2 (67%) 1 (33%)
3 (100%)
    Switch extern 230 (42%) 314 (58%)
544 (100%)
    Switch intern 37 (51%) 36 (49%)
73 (100%)
    Tussenjaar 91 (42%) 127 (58%)
218 (100%)
APCG

<0,001
    Ja 305 (46%) 364 (54%)
669 (100%)
    Nee 342 (35%) 629 (65%)
971 (100%)
    Onbekend 80 (45%) 98 (55%)
178 (100%)
Cijfer_CE_Engels 6,84 (1,13) 6,47 (1,22) <0,001 6,60 (1,21)
Cijfer_CE_Engels_missing

<0,001
    Ja 438 (45%) 539 (55%)
977 (100%)
    Nee 289 (34%) 552 (66%)
841 (100%)
Cijfer_CE_Natuurkunde 6,04 (1,00) 6,21 (0,97) 0,23 6,15 (0,98)
Cijfer_CE_Natuurkunde_missing

0,016
    Ja 637 (41%) 911 (59%)
1.548 (100%)
    Nee 90 (33%) 180 (67%)
270 (100%)
Cijfer_CE_Nederlands 6,12 (0,92) 6,17 (0,85) 0,45 6,15 (0,88)
Cijfer_CE_Nederlands_missing

<0,001
    Ja 437 (45%) 539 (55%)
976 (100%)
    Nee 290 (34%) 552 (66%)
842 (100%)
Cijfer_CE_VO 6,38 (0,37) 6,52 (0,43) <0,001 6,48 (0,41)
Cijfer_CE_VO_missing

<0,001
    Ja 414 (48%) 455 (52%)
869 (100%)
    Nee 313 (33%) 636 (67%)
949 (100%)
Cijfer_CE_Wiskunde 6,47 (1,10) 6,49 (1,18) 0,95 6,48 (1,16)
Cijfer_CE_Wiskunde_missing

<0,001
    Ja 466 (45%) 578 (55%)
1.044 (100%)
    Nee 261 (34%) 513 (66%)
774 (100%)
Cijfer_SE_VO 6,39 (0,36) 6,52 (0,43) <0,001 6,48 (0,41)
Cijfer_SE_VO_missing

<0,001
    Ja 430 (46%) 502 (54%)
932 (100%)
    Nee 297 (34%) 589 (66%)
886 (100%)
Dubbele_studie

0,066
    Ja 13 (59%) 9 (41%)
22 (100%)
    Nee 714 (40%) 1.082 (60%)
1.796 (100%)
Geslacht

<0,001
    M 177 (56%) 140 (44%)
317 (100%)
    V 550 (37%) 951 (63%)
1.501 (100%)
Leeftijd 20,75 (3,43) 19,95 (3,08) <0,001 20,27 (3,25)
Reistijd 32,18 (18,05) 32,07 (16,62) 0,50 32,11 (17,20)
SES_Arbeid -0,03 (0,10) -0,01 (0,10) <0,001 -0,02 (0,10)
SES_Totaal -0,14 (0,35) -0,07 (0,34) <0,001 -0,10 (0,35)
SES_Welvaart -0,06 (0,16) -0,03 (0,16) <0,001 -0,04 (0,16)
Studiekeuzeprofiel

<0,001
    EM 51 (40%) 78 (60%)
129 (100%)
    CM 54 (38%) 90 (63%)
144 (100%)
    EM&CM 27 (25%) 81 (75%)
108 (100%)
    NT 12 (80%) 3 (20%)
15 (100%)
    NG 175 (35%) 318 (65%)
493 (100%)
    NT&NG 50 (29%) 120 (71%)
170 (100%)
    OS 0 (0%) 3 (100%)
3 (100%)
    CERT 0 (0%) 2 (100%)
2 (100%)
    ALG 3 (50%) 3 (50%)
6 (100%)
    BI 1 (100%) 0 (0%)
1 (100%)
    EA 61 (61%) 39 (39%)
100 (100%)
    HO 9 (38%) 15 (63%)
24 (100%)
    HB 6 (55%) 5 (45%)
11 (100%)
    ICT 0 (0%) 3 (100%)
3 (100%)
    MedV 8 (47%) 9 (53%)
17 (100%)
    TP 1 (20%) 4 (80%)
5 (100%)
    TR 6 (40%) 9 (60%)
15 (100%)
    TSL 2 (67%) 1 (33%)
3 (100%)
    UV 6 (50%) 6 (50%)
12 (100%)
    VS 2 (33%) 4 (67%)
6 (100%)
    VNL 6 (60%) 4 (40%)
10 (100%)
    ZW 173 (52%) 160 (48%)
333 (100%)
Vooropleiding

<0,001
    MBO 286 (52%) 262 (48%)
548 (100%)
    HAVO 344 (35%) 636 (65%)
980 (100%)
    VWO 25 (30%) 59 (70%)
84 (100%)
    BD 47 (46%) 56 (54%)
103 (100%)
    HO 19 (26%) 53 (74%)
72 (100%)
    CD 6 (19%) 25 (81%)
31 (100%)
1 Mean (SD); n (%)
2 Wilcoxon rank sum test; Fisher’s Exact Test for Count Data with simulated p-value (based on 2000 replicates); Pearson’s Chi-squared test
Toon code 
## Toon een samenvatting van de data, gesorteerd op missende waarden
diagnose(dfOpleiding_inschrijvingen) |> 
  mutate(missing_percent = round(missing_percent, 2),
         unique_rate = round(missing_percent, 2)) |>
  arrange(desc(missing_percent)) |>
  knitr::kable(caption = "Kwaliteit van de data voor bewerkingen (gesorteerd op missende waarden)")
Kwaliteit van de data voor bewerkingen (gesorteerd op missende waarden)
variables types missing_count missing_percent unique_count unique_rate
Cijfer_CE_Natuurkunde numeric 1548 85.15 47 85.15
Cijfer_CE_Wiskunde numeric 1044 57.43 61 57.43
Cijfer_CE_Engels numeric 977 53.74 66 53.74
Cijfer_CE_Nederlands numeric 976 53.69 52 53.69
Cijfer_SE_VO numeric 932 51.27 27 51.27
Cijfer_CE_VO numeric 869 47.80 27 47.80
Studiekeuzeprofiel factor 208 11.44 23 11.44
SES_Arbeid numeric 178 9.79 311 9.79
SES_Totaal numeric 178 9.79 592 9.79
SES_Welvaart numeric 178 9.79 421 9.79
Reistijd numeric 29 1.60 379 1.60
Aanmelding numeric 0 0.00 319 0.00
Aansluiting factor 0 0.00 7 0.00
APCG character 0 0.00 3 0.00
Cijfer_CE_Engels_missing character 0 0.00 2 0.00
Cijfer_CE_Natuurkunde_missing character 0 0.00 2 0.00
Cijfer_CE_Nederlands_missing character 0 0.00 2 0.00
Cijfer_CE_VO_missing character 0 0.00 2 0.00
Cijfer_CE_Wiskunde_missing character 0 0.00 2 0.00
Cijfer_SE_VO_missing character 0 0.00 2 0.00
Collegejaar numeric 0 0.00 11 0.00
Dubbele_studie character 0 0.00 2 0.00
Geslacht character 0 0.00 2 0.00
ID character 0 0.00 1818 0.00
Leeftijd integer 0 0.00 23 0.00
Uitval factor 0 0.00 2 0.00
Vooropleiding factor 0 0.00 6 0.00

2.3 Bewerk de data

  • Uit de eerste diagnose blijkt dat niet alle variabelen goed genoeg zijn voor het bouwen van een prognosemodel: er zijn missende waarden en niet alle veldtypes zijn geschikt. We passen de variabelen aan zodat we in het model er goed mee kunnen werken.
  • Prognosemodellen kunnen niet omgaan met missende waarden. Om bias te voorkomen verwijderen we geen rijen met missende waarden, maar vullen die op (imputatie). We bewerken de data zo dat alle missende waarden worden opgevuld: bij numerieke waarden met het gemiddelde en bij categorische variabelen met ‘Onbekend’.
  • We passen sommige variabelen aan, zodat ze in het model gebruikt kunnen worden: tekstvelden zetten we om naar factor; logische variabelen (Ja/Nee) zetten we om naar een numerieke variabele (1/0).
  • De uitkomstvariabele, Uitval, leiden we af van de variabele SUC_Uitval_aantal_jaar_LTA. Als de waarde daar 1 is, is de student na 1 jaar uitgevallen, 2 na 2 jaar, etc.
  • Een fictief studentnummer (INS_Student_UUID_opleiding_vorm) gebruiken we ook, zodat we - als er afwijkende resultaten zijn - de dataset gericht kunnen onderzoeken indien nodig.
Toon code 
## Bewerk de data
dfOpleiding_inschrijvingen <- dfOpleiding_inschrijvingen |> 
  
  ## Imputeer alle numerieke variabelen met de mean
  mutate(across(where(is.numeric), ~ ifelse(
    is.na(.x),
    mean(.x, na.rm = T),
    .x
  )) ) |>
  
  ## Zet character variabelen om naar factor
  mutate(across(where(is.character), as.factor)) |> 
  
  ## Zet logische variabelen om naar 0 of 1
  mutate(across(where(is.logical), as.integer)) |>
  
  ## Vul in factoren missende waarden op met "Onbekend"
  mutate(across(where(is.factor), ~ suppressWarnings(
    fct_explicit_na(.x, na_level = "Onbekend")
  ))) |> 
  
  ## Herschik de kolommen, zodat Uitval vooraan staat
  select(Uitval, everything()) 

## Bekijk de data
## glimpse(dfOpleiding_inschrijvingen) 

## Maak een diagnose van de data
diagnose(dfOpleiding_inschrijvingen) |> 
  mutate(missing_percent = round(missing_percent, 2),
         unique_rate = round(unique_rate, 2)) |>
  knitr::kable(caption = "Kwaliteit van de data na bewerkingen")
Kwaliteit van de data na bewerkingen
variables types missing_count missing_percent unique_count unique_rate
Uitval factor 0 0 2 0.00
Aanmelding numeric 0 0 319 0.18
Aansluiting factor 0 0 7 0.00
APCG factor 0 0 3 0.00
Cijfer_CE_Engels numeric 0 0 66 0.04
Cijfer_CE_Engels_missing factor 0 0 2 0.00
Cijfer_CE_Natuurkunde numeric 0 0 47 0.03
Cijfer_CE_Natuurkunde_missing factor 0 0 2 0.00
Cijfer_CE_Nederlands numeric 0 0 52 0.03
Cijfer_CE_Nederlands_missing factor 0 0 2 0.00
Cijfer_CE_VO numeric 0 0 27 0.01
Cijfer_CE_VO_missing factor 0 0 2 0.00
Cijfer_CE_Wiskunde numeric 0 0 61 0.03
Cijfer_CE_Wiskunde_missing factor 0 0 2 0.00
Cijfer_SE_VO numeric 0 0 27 0.01
Cijfer_SE_VO_missing factor 0 0 2 0.00
Collegejaar numeric 0 0 11 0.01
Dubbele_studie factor 0 0 2 0.00
Geslacht factor 0 0 2 0.00
ID factor 0 0 1818 1.00
Leeftijd integer 0 0 23 0.01
Reistijd numeric 0 0 379 0.21
SES_Arbeid numeric 0 0 311 0.17
SES_Totaal numeric 0 0 592 0.33
SES_Welvaart numeric 0 0 421 0.23
Studiekeuzeprofiel factor 0 0 23 0.01
Vooropleiding factor 0 0 6 0.00

2.4 Bekijk de onderlinge correlaties

Het is verstandig om voorafgaand aan het bouwen van een model te kijken naar de onderlinge correlaties tussen numerieke variabelen. Dit geeft inzicht in de data en kan helpen bij het maken van keuzes voor het model of de duiding van de uitkomsten.

Toon code 
## Maak een plot van de onderlinge correlaties in numerieke variabelen
dfOpleiding_inschrijvingen |> 
  select(-Collegejaar) |>
  select(where(is.numeric)) |> 
  cor() |> 
  corrplot::corrplot(
    order = 'hclust', 
    addrect = 4,
    method = "number",  
    tl.cex = 0.8,       
    tl.col = "black",
    diag = FALSE)

2.5 Bouw de trainingset, testset en validatieset

  • De data is nu geschikt om een prognosemodel mee te bouwen.
  • Om het model te bouwen, testen en valideren, splitsen we de data in drie delen van 60%, 20% en 20%. We doen dit op zo’n manier, dat elk deel ongeveer een gelijk aantal studenten bevat dat uitvalt.
  • We trainen het model op basis van 60% en testen de modellen tijden het trainen op de overige 20% (de testset).
  • Als het model klaar is, valideren we het op de 20% studenten uit de validatieset, die is opgebouwd uit zo recent mogelijke data. De validatieset blijft dus de gehele tijd ongemoeid, zodat we overfitting - een te goed model op bekende data, maar slechte performance op onbekende data - voorkomen.
  • Een willekeurig, maar vaststaand seedgetal voorkomt dat we bij elke run van het model c.q. deze code een net iets andere uitkomst krijgen.
Toon code 
set.seed(0821)

## Splits de data in 3 delen: 60%, 20% en 20%
splits      <- initial_validation_split(dfOpleiding_inschrijvingen,
                                        strata = Uitval,
                                        prop = c(0.6, 0.2))

## Maak drie sets: een trainingset, een testset en een validatieset
dfUitval_train      <- training(splits)
dfUitval_test       <- testing(splits)
dfUitval_validation <- validation_set(splits)

## Maak een resample set op basis van 10 folds (default)
dfUitval_resamples  <- vfold_cv(dfUitval_train, strata = Uitval)

## Training set proporties uitval
dfUitval_train |> 
  count(Uitval) |> 
  mutate(prop = n/sum(n)) |> 
  knitr::kable(caption = "Trainingsset")
Trainingsset
Uitval n prop
FALSE 654 0.6
TRUE 436 0.4
Toon code 
## Test set proporties uitval
dfUitval_test  |> 
  count(Uitval) |> 
  mutate(prop = n/sum(n)) |> 
  knitr::kable(caption = "Testset")
Testset
Uitval n prop
FALSE 219 0.6
TRUE 146 0.4

3 Model I: Logistische Regressie

  • Het eerste model is een logistische regressie met penalized likelihood; we gebruiken de glmnet engine voor het bouwen van het model. Penalized likelihood is een techniek die helpt bij het voorkomen van overfitting. Glmnet is een populair package voor het bouwen van logistische regressiemodellen.
  • We gebruiken de Area under the ROC Curve (AUC/ROC) als performance metric.

3.1 Maak het model

We bouwen eerst het model.

## Bouw het model: logistische regressie
lr_mod <- 
  logistic_reg(penalty = tune(), mixture = 1) |> 
  set_engine("glmnet")

3.2 Maak de recipe

Vervolgens zetten we meerdere stappen in een ‘recipe’. We verwijderen de student ID en het collegejaar uit de data, omdat deze niet moet worden gebruikt in het model. We converteren factoren naar dummy variabelen, verwijderen zero values en centreren en schalen numerieke variabelen.

## Bouw de recipe: logistische regressie
lr_recipe <- 
  recipe(Uitval ~ ., data = dfUitval_train) |>  
  update_role(ID, new_role = "ID") |>           ## Zet de student ID als ID variabele
  step_rm(ID, Collegejaar) |>                   ## Verwijder ID en collegejaar uit het model
  step_dummy(all_nominal_predictors()) |>       ## Converteer factoren naar dummy variabelen
  step_zv(all_predictors()) |>                  ## Verwijder zero values
  step_normalize(all_numeric_predictors())      ## Centreer en schaal numerieke variabelen

## Toon de recipe
tidy(lr_recipe) |> 
  knitr::kable()
number operation type trained skip id
1 step rm FALSE FALSE rm_yTMrb
2 step dummy FALSE FALSE dummy_0LqEE
3 step zv FALSE FALSE zv_P3EsP
4 step normalize FALSE FALSE normalize_Yy34e

3.3 Maak de workflow

Voor de uitvoering bouwen we een nieuwe workflow. Daaraan voegen we het model en de bewerkingen in de recipe toe.

## Maak de workflow: logistische regressie
lr_workflow <- 
  workflow() |>         ## Maak een workflow
  add_model(lr_mod) |>  ## Voeg het model toe
  add_recipe(lr_recipe) ## Voeg de recipe toe

## Toon de workflow
lr_workflow
══ Workflow ════════════════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: logistic_reg()

── Preprocessor ────────────────────────────────────────────────────────────────
4 Recipe Steps

• step_rm()
• step_dummy()
• step_zv()
• step_normalize()

── Model ───────────────────────────────────────────────────────────────────────
Logistic Regression Model Specification (classification)

Main Arguments:
  penalty = tune()
  mixture = 1

Computational engine: glmnet

3.4 Tune en train het model

Het model moet getuned worden. Dit houdt in dat we de beste parameters voor het model moeten vinden. We maken een grid met verschillende penalty waarden. Daarmee kunnen we vervolgens het beste model selecteren met de hoogste ROC/AUC. We plotten de resultaten van de tuning, zodat we hieruit het beste model kunnen kiezen.

## Maak een grid: logistische regressie
lr_reg_grid <- tibble(penalty = 10 ^ seq(-4, -1, length.out = 30))

## Train en tune het model: logistische regressie
lr_res <- 
  lr_workflow |> 
  tune_grid(dfUitval_validation,
            grid = lr_reg_grid,
            control = control_grid(save_pred = TRUE),
            metrics = metric_set(roc_auc))
Toon code 
## Plot de resultaten + een rode verticale lijn voor de max AUC
lr_plot <- 
  lr_res |> 
  collect_metrics() |> 
  ggplot(aes(x = penalty, y = mean)) + 
  geom_point() + 
  geom_line() + 
  ylab("Area under the ROC Curve") +
  scale_x_log10(labels = scales::label_number())

# Zoek de penalty waarde met de max AUC
max_auc_penalty <- lr_res |> 
  collect_metrics() |> 
  filter(mean == max(mean)) |> 
  pull(penalty)

# Voeg de rode verticale lijn toe aan lr_plot
lr_plot_plus <- lr_plot + 
  geom_vline(xintercept = max_auc_penalty, color = "red")

# Vindt een mean voor de max AUC die hoger is
max_auc_mean <- lr_res |> 
  collect_metrics() |> 
  filter(mean == max(mean)) |> 
  pull(penalty)


lr_plot_plus

3.5 Kies het beste model

We evalueren modellen met een zo hoog mogelijke Area under the ROC Curve (AUC/ROC) en een zo laag mogelijke penalty. Zo kunnen we uit de resultaten het beste model kiezen. Tot slot maken we een ROC curve om de prestaties van het model te visualiseren.

## Toon het beste model
top_models <-
  lr_res |> 
  show_best(metric = "roc_auc", n = 10) |> 
  mutate(mean = round(mean, 6)) |>
  arrange(penalty) 

top_models|> 
  knitr::kable()
penalty .metric .estimator mean n std_err .config
0.0035622 roc_auc binary 0.652199 1 NA Preprocessor1_Model16
0.0045204 roc_auc binary 0.652958 1 NA Preprocessor1_Model17
0.0057362 roc_auc binary 0.654160 1 NA Preprocessor1_Model18
0.0072790 roc_auc binary 0.654919 1 NA Preprocessor1_Model19
0.0092367 roc_auc binary 0.655489 1 NA Preprocessor1_Model20
0.0117210 roc_auc binary 0.655394 1 NA Preprocessor1_Model21
0.0148735 roc_auc binary 0.653243 1 NA Preprocessor1_Model22
0.0188739 roc_auc binary 0.655220 1 NA Preprocessor1_Model23
0.0239503 roc_auc binary 0.656675 1 NA Preprocessor1_Model24
0.0303920 roc_auc binary 0.655315 1 NA Preprocessor1_Model25 
## Selecteer het beste model: logistische regressie
lr_best <- 
  lr_res |> 
  collect_metrics() |> 
  filter(mean == max(mean)) |>
  slice(1) 

lr_best|> 
  mutate(mean = round(mean, 6)) |>
  knitr::kable()
penalty .metric .estimator mean n std_err .config
0.0239503 roc_auc binary 0.656675 1 NA Preprocessor1_Model24 
## Verzamel de predicties en evalueer het model (AUC/ROC): logistische regressie
lr_auc <- 
  lr_res |> 
  collect_predictions(parameters = lr_best) |> 
  roc_curve(Uitval, .pred_FALSE) |> 
  mutate(model = "Logistisch Regressie")

autoplot(lr_auc) 

## Bepaal de AUC van het beste model
lr_auc_highest   <-
  lr_res |>
  collect_predictions(parameters = lr_best) |> 
  roc_auc(Uitval, .pred_FALSE)

## Voeg de naam van het model en de AUC toe dfModel_results
dfModel_results <- 
  dfModel_results |>
  add_row(model = "Logistic Regression", auc = lr_auc_highest$.estimate)

4 Model II: Tree-based ensemble

  • Het tweede model is een random forest: een ensemble van beslisbomen (decision trees). Het is een krachtig model dat goed om kan gaan met complexe data en veel variabelen.
  • We gebruiken de ranger engine voor het bouwen van het model.

4.1 Bepaal het aantal PC-cores

Omdat een random forest model veel berekeningen vereist, willen we daarvoor alle computerkracht gebruiken die beschikbaar is. Het aantal CPU’s (cores) van de computer bepaalt hoe snel het model getraind kan worden. Deze informatie gebruiken we bij het bouwen van het model.

Toon code 
## Bepaal het aantal cores
cores <- parallel::detectCores()

4.2 Maak het model

We bouwen eerst het model. We gebruiken de rand_forest functie om het model te bouwen. We tunen de mtry en min_n parameters. De mtry parameter bepaalt het aantal variabelen dat per boom wordt gebruikt. De min_n parameter bepaalt het minimum aantal observaties dat in een blad van de boom moet zitten. De functie tune() is hier nog een placeholder om de beste waarden voor deze parameters - die we later bepalen - daar in te stellen. We gebruiken 1.000 bomen c.q. versies van het model.

## Bouw het model: random forest

rf_mod <- 
  rand_forest(mtry = tune(), min_n = tune(), trees = 1000) |> 
  set_engine("ranger", num.threads = cores) |> 
  set_mode("classification")

4.3 Maak de recipe

We maken een recipe voor het random forest model. We verwijderen de student ID en het collegejaar uit de data, omdat deze niet moet worden gebruikt in het model. Overige stappen zijn bij een random forest minder relevant in tegenstelling tot een regressiemodel.

## Maak de recipe: random forest
rf_recipe <- 
  recipe(Uitval ~ ., data = dfUitval_train) |> 
  step_rm(ID, Collegejaar)                      ## Verwijder ID en Collegejaar uit het model
  
## Toon de recipe
tidy(rf_recipe) |> 
  knitr::kable()
number operation type trained skip id
1 step rm FALSE FALSE rm_u6aow

4.4 Maak de workflow

We voegen het model en de recipe toe aan de workflow voor dit model.

## Maak de workflow: random forest
rf_workflow <- 
  workflow() |> 
  add_model(rf_mod) |> 
  add_recipe(rf_recipe)

## Toon de workflow
rf_workflow
══ Workflow ════════════════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: rand_forest()

── Preprocessor ────────────────────────────────────────────────────────────────
1 Recipe Step

• step_rm()

── Model ───────────────────────────────────────────────────────────────────────
Random Forest Model Specification (classification)

Main Arguments:
  mtry = tune()
  trees = 1000
  min_n = tune()

Engine-Specific Arguments:
  num.threads = cores

Computational engine: ranger

4.5 Tune en train het model

We trainen en tunen het model in de workflow. We maken een grid met verschillende waarden voor de parameters mtry en min_n. We gebruiken de Area under the ROC Curve (AUC/ROC) als performance metric. Met de resultaten van de tuning kiezen we het beste model.

## Toon de parameters die getuned kunnen worden
rf_mod
Random Forest Model Specification (classification)

Main Arguments:
  mtry = tune()
  trees = 1000
  min_n = tune()

Engine-Specific Arguments:
  num.threads = cores

Computational engine: ranger 
## Extraheer de parameters die getuned worden
extract_parameter_set_dials(rf_mod)
Collection of 2 parameters for tuning

 identifier  type    object
       mtry  mtry nparam[?]
      min_n min_n nparam[+]

Model parameters needing finalization:
   # Randomly Selected Predictors ('mtry')

See `?dials::finalize` or `?dials::update.parameters` for more information.
## Bepaal de seed
set.seed(2904)

## Bouw het grid: random forest
rf_res <- 
  rf_workflow |> 
  tune_grid(dfUitval_validation,
            grid = 25,
            control = control_grid(save_pred = TRUE),
            metrics = metric_set(roc_auc))
i Creating pre-processing data to finalize unknown parameter: mtry

4.6 Kies het beste model

We evalueren de beste modellen en maken een ROC curve om de performance van het model te visualiseren. Vervolgens vergelijken we de prestaties van de modellen en kiezen we het beste model.

## Toon de beste modellen
rf_res |> 
  show_best(metric = "roc_auc", n = 15) |> 
  mutate(mean = round(mean, 6)) |>
  knitr::kable()
mtry min_n .metric .estimator mean n std_err .config
9 12 roc_auc binary 0.655394 1 NA Preprocessor1_Model05
18 40 roc_auc binary 0.653464 1 NA Preprocessor1_Model16
11 29 roc_auc binary 0.651693 1 NA Preprocessor1_Model02
2 6 roc_auc binary 0.651534 1 NA Preprocessor1_Model22
5 36 roc_auc binary 0.650870 1 NA Preprocessor1_Model09
16 34 roc_auc binary 0.650617 1 NA Preprocessor1_Model07
4 26 roc_auc binary 0.650585 1 NA Preprocessor1_Model23
13 25 roc_auc binary 0.650554 1 NA Preprocessor1_Model25
3 8 roc_auc binary 0.650364 1 NA Preprocessor1_Model08
7 35 roc_auc binary 0.650142 1 NA Preprocessor1_Model14
20 10 roc_auc binary 0.650142 1 NA Preprocessor1_Model18
11 21 roc_auc binary 0.649953 1 NA Preprocessor1_Model15
15 20 roc_auc binary 0.649541 1 NA Preprocessor1_Model20
22 38 roc_auc binary 0.649415 1 NA Preprocessor1_Model01
17 24 roc_auc binary 0.649320 1 NA Preprocessor1_Model10 
## Plot de resultaten
autoplot(rf_res) 

## Selecteer het beste model
rf_best <- 
  rf_res |> 
  select_best(metric = "roc_auc")

rf_best|> 
  knitr::kable()
mtry min_n .config
9 12 Preprocessor1_Model05
Toon code 
## Verzamel de predicties
rf_res |> 
  collect_predictions() |> 
  head(10) |>
  knitr::kable()
.pred_FALSE .pred_TRUE id .row mtry min_n Uitval .config
0.6355015 0.3644985 validation 1091 22 38 TRUE Preprocessor1_Model01
0.7976896 0.2023104 validation 1092 22 38 FALSE Preprocessor1_Model01
0.5274335 0.4725665 validation 1093 22 38 FALSE Preprocessor1_Model01
0.8106095 0.1893905 validation 1094 22 38 FALSE Preprocessor1_Model01
0.6987968 0.3012032 validation 1095 22 38 FALSE Preprocessor1_Model01
0.3413777 0.6586223 validation 1096 22 38 TRUE Preprocessor1_Model01
0.9585159 0.0414841 validation 1097 22 38 FALSE Preprocessor1_Model01
0.5565431 0.4434569 validation 1098 22 38 TRUE Preprocessor1_Model01
0.8452084 0.1547916 validation 1099 22 38 TRUE Preprocessor1_Model01
0.7423587 0.2576413 validation 1100 22 38 FALSE Preprocessor1_Model01
Toon code 
## Bepaal de AUC/ROC curve
rf_auc <- 
  rf_res |> 
  collect_predictions(parameters = rf_best) |> 
  roc_curve(Uitval, .pred_FALSE) |> 
  mutate(model = "Random Forest")

autoplot(rf_auc)

Toon code 
## Bepaal de AUC van het beste model
rf_auc_highest   <-
  rf_res |>
  collect_predictions(parameters = rf_best) |> 
  roc_auc(Uitval, .pred_FALSE)

## Voeg de naam van het model en de AUC toe dfModel_results
dfModel_results <- 
  dfModel_results |>
  add_row(model = "Random Forest", auc = rf_auc_highest$.estimate)

5 De uiteindelijke fit

  • In de laatste stap van deze analyse maken we het model definitief.
  • We testen het model op de testset en evalueren het model met metrieken en de Variable Importance Factor (VIF).

5.1 Combineer de AUC/ROC curves en kies het beste model

Eerst combineren we de AUC/ROC curves van de modellen om ze te vergelijken. We kiezen het beste model op basis van de hoogste AUC/ROC.

Toon code 
## Combineer de AUC/ROC curves om de modellen te vergelijken
bind_rows(lr_auc, rf_auc) |> 
  ggplot(aes(x = 1 - specificity, y = sensitivity, col = model)) + 
  geom_path(lwd = 1.5, alpha = 0.8) +
  geom_abline(lty = 3) + 
  coord_equal() + 
  scale_color_viridis_d(option = "plasma", end = .6)

Toon code 
## Bepaal welke van de modellen het beste is op basis van de hoogste AUC/ROC
dfModel_results <- dfModel_results |>
  mutate(number = row_number()) |> 
  mutate(best = ifelse(auc == max(auc), TRUE, FALSE)) |> 
  arrange(number)

## Bepaal het beste model
sBest_model <- dfModel_results$model[dfModel_results$best == TRUE]
sBest_model_auc <- round(dfModel_results$auc[dfModel_results$best == TRUE], 4)

Het beste model is het Logistic Regression model met een AUC/ROC van 0.6567. Het Logistic Regression model heeft een AUC van 0.6567. Het Random Forest model heeft een AUC van 0.6554. We ronden de analyse verder af met het Logistic Regression model.

5.2 Maak het finale model

We maken het finale model op basis van de beste parameters die we hebben gevonden. Door in de engine bij importance de impurity op te geven, wordt het beste random forest model gekozen om de data definitief mee te classificeren.

## Test het ontwikkelde model op de testset
## Bepaal de optimale parameters

## Bouw het laatste model
if (sBest_model == "Logistic Regression") {
  
  last_lr_mod <-
    logistic_reg(penalty = lr_best$penalty,
                 mixture = 1) |>
    set_engine("glmnet") |>
    set_mode("classification")

} else if (sBest_model == "Random Forest") {
  
  last_rf_mod <-
    rand_forest(mtry = rf_best$mtry,
                min_n = rf_best$min_n,
                trees = 1000) |>
    set_engine("ranger", num.threads = cores, importance = "impurity") |>
    set_mode("classification")

}

5.3 Maak de workflow

We voegen het model toe aan de workflow en updaten de workflow met het finale model.

## Bouw de laatste workflow op basis van het laatste model
if (sBest_model == "Logistic Regression") {
  
  last_lr_workflow <- 
    lr_workflow |> 
    update_model(last_lr_mod)
  
} else if (sBest_model == "Random Forest") {
  
  last_rf_workflow <- 
    rf_workflow |> 
    update_model(last_rf_mod)
  
}

5.4 Fit het finale model

We voeren de finale fit uit. De functie last_fit past het model toe op de validatieset.

## Voer de laatste fit uit
set.seed(2904)

if(sBest_model == "Logistic Regression") {
  
  last_fit <- 
    last_lr_workflow |> 
    last_fit(splits)
  
} else if(sBest_model == "Random Forest") {
  
  last_fit <- 
    last_rf_workflow |> 
    last_fit(splits)
  
}

5.5 Evalueer het finale model: metrieken en vif

We evalueren het finale model op basis van 3 metrieken (accuraatheid, ROC/AUC en de Brier score = de Mean Squared Error) en de Variable Importance Factor (VIF). Uit de VIF is op te maken welke variabelen het meest bijdragen aan de voorspelling van de uitkomstvariabele.

## Verzamel de metrieken
last_fit |> 
  collect_metrics() |> 
  mutate(.estimate = round(.estimate, 4)) |>
  knitr::kable()
.metric .estimator .estimate .config
accuracy binary 0.6548 Preprocessor1_Model1
roc_auc binary 0.6457 Preprocessor1_Model1
brier_class binary 0.2246 Preprocessor1_Model1 
## Extraheer de feature importance
last_fit |>
  extract_fit_parsnip() |>
  vip(num_features = 20) +
  ggtitle("Meest voorspellende factoren") +
  ylab("Importance")

6 Verdiepende analyse van het model

We weten nu welke variabelen van invloed zijn, maar niet hoe en in welke richting: dragen ze sterk bij of juist niet, verhogen of verlagen ze uitval? Om het model beter te begrijpen en te kunnen uitleggen, maken met behulp van het Dalex package een explainer. Dit package is ontwikkeld om beter uit te leggen welke variabelen van belang zijn en wat deze voor een effect hebben in een model. Een explainer is een model-onafhankelijke wrapper om het model heen en geeft inzicht in de voorspellingen van het model en de bijdrage van de variabelen aan de prognose. Een explainer maakt het verder mogelijk om modellen onderling te vergelijken en benchmarken.

6.1 Maak een explainer

Op basis van het tidymodels model extraheren we de informatie voor de explainer.

Toon code 
fitted_model_tmp <- last_fit |>
  extract_fit_parsnip()

vi_df <- vip::vi(fitted_model_tmp)
vi_df |> 
  head(20) |> 
  knitr::kable()
Variable Importance Sign
Cijfer_CE_VO 1.0779172 NEG
Aansluiting_Direct 0.6245274 POS
Cijfer_CE_VO_missing_Nee 0.5570801 NEG
Aansluiting_Switch.extern 0.5480475 POS
Cijfer_SE_VO 0.5100031 POS
Cijfer_SE_VO_missing_Nee 0.4853102 POS
Aansluiting_Tussenjaar 0.4456396 POS
APCG_Onbekend 0.3628895 POS
Cijfer_CE_Nederlands_missing_Nee 0.3591258 POS
Cijfer_CE_Engels 0.3429742 POS
Geslacht_V 0.3208378 NEG
Aanmelding 0.3156921 NEG
Cijfer_CE_Engels_missing_Nee 0.2954474 NEG
Studiekeuzeprofiel_OS 0.2894741 NEG
Studiekeuzeprofiel_ZW 0.2578913 POS
Studiekeuzeprofiel_EM.CM 0.2524596 NEG
SES_Welvaart 0.2468578 NEG
Studiekeuzeprofiel_ICT 0.2344931 NEG
Studiekeuzeprofiel_NG 0.2284218 NEG
Aansluiting_Switch.intern 0.2257438 POS
Toon code 
## Extraheer het fitted model en de workflow
fitted_model <- last_fit |>
  extract_fit_parsnip()

workflow <- last_fit |>
  extract_workflow()

# Pas de Uitval variabele aan naar numeric (0/1), 
# zodat er een explainer van gemaakt kan worden
dfOpleiding_inschrijvingen$Uitval <- as.numeric(dfOpleiding_inschrijvingen$Uitval) - 1

# Maak een explainer
explain_lm <- DALEX::explain(
  model = workflow,
  data = dfOpleiding_inschrijvingen,
  y = dfOpleiding_inschrijvingen$Uitval,
  label = "Linear Regression")
Preparation of a new explainer is initiated
  -> model label       :  Linear Regression 
  -> data              :  1818  rows  27  cols 
  -> target variable   :  1818  values 
  -> predict function  :  yhat.workflow  will be used (  default  )
  -> predicted values  :  No value for predict function target column. (  default  )
  -> model_info        :  package tidymodels , ver. 1.2.0 , task classification (  default  ) 
  -> predicted values  :  numerical, min =  0.09233298 , mean =  0.3963696 , max =  0.738768  
  -> residual function :  difference between y and yhat (  default  )
  -> residuals         :  numerical, min =  -0.7024894 , mean =  0.003520435 , max =  0.8469217  
  A new explainer has been created!
Toon code 
library(iBreakDown)

## Bereken de model parts op basis van de RMSE
mp <- model_parts(explain_lm, loss_function = loss_root_mean_square)

# Toon de resultaten
plot(mp)

6.2 Pas het model toe op een willekeurige student

Nu we deze explainer hebben, kunnen we het model toepassen op een willekeurige student. We kiezen een student uit de dataset en passen het model toe op deze student. We kijken naar de voorspelling van het model en de bijdrage van de variabelen aan deze voorspelling. Dit geeft een verder inzicht in de werking van het model.

De intercept van het model (de basisvoorspelling) is 0,396: een kans op uitval van 39.6%. De bijdrage van de variabelen aan de voorspelling kan positief of negatief zijn. Een positieve bijdrage betekent dat de variabele de voorspelling verhoogt, een negatieve bijdrage betekent dat de voorspelling wordt verlaagd. De combinatie van variabelen bij deze student verlaagt de kans op uitval naar 0,254 (25,4%).

Toon code 
source("01. Includes/Studentpersonas_opleidingen.R")

## Loop over de persona's
for (i in 1:nrow(dfPersona)) {

  student_x <- dfPersona[i,]

  bd_lm <- predict_parts(explainer = explain_lm,
                     new_observation = student_x,
                                type = "break_down")

  # Convert to data frame for ggplot
  bd_lm_df <- as.data.frame(bd_lm) |> 

    ## Sorteer op basis van de position
    arrange(desc(position)) |> 

    ## Maak het label aan
    mutate(label = paste0(as.character(round(contribution, 3)*100),"%")) |> 
    mutate(label = case_when(
      variable %in% c("intercept","prediction") ~ label,
      sign == 1  ~ paste0("+", label),
      sign == -1 ~ paste0(label),
      sign == 0  ~ paste0("+", label),
      .default = label
    )) |> 
    mutate(label = case_when(
      label == "0%" ~ "+0%",
      .default = label
    )) |>

  filter(sign != 0) |>

  mutate(start = dplyr::lag(cumulative, default = min(cumulative)),
         end = cumulative) 

  ## Voeg een rij toe "Overige variabelen"
  bd_lm_df <- bd_lm_df |> 
     add_row(variable = "+ Overige variabelen", 
          contribution = 0, 
          variable_name = NA,
          variable_value = NA,
          cumulative = bd_lm_df$cumulative[bd_lm_df$position == 1],
          sign = "0", 
          position = 2,
          label = "+0%", 
          start = bd_lm_df$cumulative[bd_lm_df$position == 1], 
          end = bd_lm_df$cumulative[bd_lm_df$position == 1]
          ) |> 
  arrange(position) |>
  mutate(position = row_number()) |> 
  arrange(desc(position)) |> 
  mutate(next_start = lead(start, default = NA),
         next_position = lead(position, default = NA)) |> 
  mutate(start = case_when(
          variable == "intercept" ~ 0,
          variable == "prediction" ~ bd_lm_df$end[bd_lm_df$variable == "intercept"],
          .default = start)) |> 
  mutate(sign = case_when(
          variable == "intercept" ~ "X",
          .default = sign)) |> 
  mutate(label_color = case_when(
          sign == "1" ~ "#C8133B",
          sign == "-1" ~ "#466F9D",
          .default = "black")) |> 
  rowwise() |>
  mutate(label_position = max(start, end)) |> 
  ungroup()

  bd_lm_df$sign <- factor(bd_lm_df$sign, 
                          levels = c("1", "-1", "0", "X"), 
                          labels = c("Positief", "Negatief", "Geen", "X")) 
  
  nUitval <- round(bd_lm_df$cumulative[bd_lm_df$variable == 'prediction'], 3)*100
  student_x_title <- paste0("Gemiddelde student | vooropleiding: ", student_x$Vooropleiding,
                            " | kans op uitval: ",
                            nUitval, "%")

  # Bepaal de breaks voor de x-as
  y_breaks <- seq(0, 1, by = 0.2)
  y_labels <- paste0(seq(0, 100, by = 20), "%")

  # Watervalplot met ggplot
  p <- ggplot(bd_lm_df) +
    # Voeg horizontale lijnen toe om de 0.2
    geom_hline(yintercept = y_breaks, color = "#CBCBCB", linetype = "solid", linewidth = 0.5) +
    geom_hline(yintercept = min(bd_lm_df$cumulative), color = "black", linetype = "dotted") +
    geom_rect(aes(xmin = position - 0.4, xmax = position + 0.4,
                  ymin = start, ymax = end, fill = sign)) +
    geom_segment(aes(x = next_position - 0.4, xend = position + 0.4,
                     y = end, yend = end), color = "darkgray") +
    coord_flip() +
    labs(title = "Breakdown van kans op Uitval",
         subtitle = student_x_title,
         x = NULL,
         y = NULL) +
    scale_fill_manual(values = c("Positief" = "#C8133B", "Negatief" = "#466F9D", "X" = "#57606C")) +  
    # Voeg tekstlabels toe voor de variabelen
    geom_text(aes(x = position, y = label_position, label = label),
              hjust = -0.1, 
              size = 4,
              color = "black") +
    
    # Pas het thema aan om de y-as labels weer te geven
    scale_x_continuous(breaks = bd_lm_df$position, labels = bd_lm_df$variable) +
    scale_y_continuous(breaks = y_breaks, labels = y_labels, limits = c(0, 1)) +
    theme_minimal() +
    theme(legend.position = "none") +
    theme(axis.text.y = element_text(size = 10),
          panel.grid.major = element_blank(),  
          panel.grid.minor = element_blank() ) 
  
  # Print de plot
  suppressWarnings(print(p))

}

6.3 Plot de ROC curve

Tot slot maken we een ROC curve om de prestaties van het definitieve model te visualiseren. De Sensitivity (True Positive Rate) en Specificity (True Negative Rate) worden hierin uitgezet. De Area under the ROC Curve (AUC/ROC) geeft de prestaties van het model weer. Het model scoort beter naarmate de AUC/ROC dichter bij de 1 ligt (in de linkerbovenhoek). Een AUC/ROC van 0,5 betekent dat het model niet beter presteert dan een willekeurige voorspelling.

## Toon de roc curve
last_fit |> 
  collect_predictions() |> 
  roc_curve(Uitval, .pred_FALSE) |> 
  autoplot()

7 Conclusies

7.1 Het beste prognosemodel voor deze opleiding

Het beste prognosemodel blijkt het Logistic Regression model te zijn.

  • Van de prognosemodellen die we hebben ontwikkeld om uitval na 1 jaar te voorspellen, had het Logistic Regression model de hoogste AUC/ROC waarde (0.6567).

7.2 Mate van accuraatheid en lift

Een prognosemodel moet minimaal beter presteren dan een base-model om waarde toe te voegen. Het base-model neemt de grootste klasse van de gemiddelde uitval na 1 jaar van de afgelopen jaren als basis. Stel we zouden tegen alle studenten zeggen dat ze hun studie gaan halen, dan is de mate van accuratesse gelijk aan dit base-model. Dit base-model is dus altijd hoger dan de 50% lijn van de AUC/ROC curve, tenzij het base-model toevallig precies 50% is

De mate van accuraatheid van de toepassing van het model is vrij laag (65.48%).

  • Base-model: 60.01% – Voor deze opleiding bereken we het base-model als volgt. Van alle studenten viel 39.99% uit. De grootste klasse (en de accuratesse) van het base-model is daarmee (100% - 39.99% = ) 60.01% die niet uitviel.
  • Accuratesse prognose: 65.48% – Het model voorspelt Uitval na 1 jaar met een accuratesse van 65.48%.
  • Lift: 5.47% – Het model scoort in de huidige opbouw met een verschil van 5.47% (de lift) wat beter dan de accuraatheid van het base-model (60.01%).

7.3 Factoren

  • Vooralsnog geen nadere bijzonderheden.

 

Verantwoording

Deze analyse maakt deel uit van het onderzoek naar kansengelijkheid van het lectoraat Learning Technology & Analytics van De Haagse Hogeschool: No Fairness without Awareness | Het rapport is door het lectoraat ontwikkeld in Quarto 1.4.549. | Template versie: 0.9.1.9000

 

Copyright

Dr. Theo Bakker, Lectoraat Learning Technology & Analytics, De Haagse Hogeschool © 2023-2024 Alle rechten voorbehouden.

Warning in rm(rf_auc, rf_auc_highest, rf_best, rf_mod, rf_recipe, rf_res, :
object 'last_rf_mod' not found
Warning in rm(rf_auc, rf_auc_highest, rf_best, rf_mod, rf_recipe, rf_res, :
object 'last_rf_workflow' not found